namespace NBody
{
    using System;
    using System.IO;
    using System.Runtime.InteropServices;
    
    public class Particle : IComparable
    {
        protected double mass;
        protected double x;
        protected double y;
        protected double z;
        protected double vx;
        protected double vy;
        protected double vz;
        
        public Particle(double mass, double x, double y, double z,
                 double vx, double vy, double vz)
        {
            this.mass = mass;
            this.x = x;
            this.y = y;
            this.z = z;
            this.vx = vx;
            this.vy = vy;
            this.vz = vz;
        }
        
        public Particle(Particle p) : this(p.Mass, p.X, p.Y, p.Z, p.Vx, p.Vy, p.Vz)
        {
        }
        
        public double Mass { get { return mass; } }
                
        public double X { get { return x; } }
        public double Y { get { return y; } }
        public double Z { get { return z; } }
        
        public double Position(int i)
        {
            if (i == 0)
                return x;
            else if (i == 1)
                return y;
            else if (i == 2)
                return z;
            else
                throw new NBodyException(String.Format("Particle.Position({0}) not valid", i));
        }
        
        public double[] PositionArray
        {
            get { return new double[] {x, y, z} ;}
        }
        
        public double Velocity(int i)
        {
            if (i == 0)
                return vx;
            else if (i == 1)
                return vy;
            else if (i == 2)
                return vz;
            else
                throw new NBodyException(String.Format("Particle.Velocity({0}) not valid", i));
        }
        
        public double[] VelocityArray
        {
            get { return new double[] {x, y, z} ;}
        }
        
        public void SetMass(double m)
        {
            this.mass = m;
        }
        
        public void SetPosition(double x, double y, double z)
        {
            this.x = x;
            this.y = y;
            this.z = z;
        }
        
        public double Vx { get { return vx; } }
        public double Vy { get { return vy; } }
        public double Vz { get { return vz; } }
        
        public void SetVelocity(double vx, double vy, double vz)
        {
            this.vx = vx;
            this.vy = vy;
            this.vz = vz;
        }
        
        public double Radius
        {
            get
            {
                return Math.Sqrt(X * X + Y * Y + Z * Z);
            }
        }
        
        public double RadiusSquared
        {
            get
            {
                return X * X + Y * Y + Z * Z;
            }
        }
        
        public double AbsoluteVelocity
        {
            get
            {
                return Math.Sqrt(Vx * Vx + Vy * Vy + Vz * Vz);
            }
        }
        
        public double VelocitySquared
        {
            get
            {
                return Vx * Vx + Vy * Vy + Vz * Vz;
            }
        }
        
        public double RadialVelocity
        {
            get
            {
                double r = Radius;
                return vx * x / r + vy * y / r + vz * z / r;
            }
        }
        
        public double ThetaVelocity
        {
            get
            {
                return vx * Math.Cos(Theta) * Math.Cos(Phi) + 
                       vy * Math.Cos(Theta) * Math.Sin(Phi) -
                       vz * Math.Sin(Theta);
            }
        }
        
        public double PhiVelocity
        {
            get
            {
                return -vx * Math.Sin(Phi) + vy * Math.Cos(Phi);
            }
        }
        
        public double TangentialVelocity
        {
            get
            {
                double r = Radius;
                double xx = y * vz - z * vy;
                double yy = z * vx - x * vz;
                double zz = x * vy - y * vx;
                return Math.Sqrt(xx*xx + yy*yy + zz*zz) / r;
            }
        }
        
        /* yuck -- although this is correct, it leads to poor thinking
           about J-vector.  Suggest using J^2, along with Jx, Jy, Jz
        public double AngularMomentum // per unit mass
        {
            get
            {
                return TangentialVelocity * Radius;
            }
        }
        */
                
        // angular momentum squared (per unit mass)
        public double J2
        {
            get
            {
                double xx = y * vz - z * vy;
                double yy = z * vx - x * vz;
                double zz = x * vy - y * vx;
                return xx * xx + yy * yy + zz * zz;
            }
        }
        
        // x-component of angular momentum
        public double Jx
        {
            get { return y * vz - z * vy; }
        }
        
        // y-component of angular momentum
        public double Jy
        {
            get { return z * vx - x * vz; }
        }
        
        // z-component of angular momentum
        public double Jz
        {
            get { return x * vy - y * vx; }
        }
        
        // spherical coord theta (angle between position vector
        // and the z-axis
        public double Theta
        {
            get
            {
                return Math.Acos(z / Radius);
            }
        }
        
        // spherical coord phi (azimuthal angle)
        public double Phi
        {
            get
            {
                return Math.Atan2(y, x);
            }
        }
        
        public double CylindricalRadius
        {
            get
            {
                return Math.Sqrt(x * x + y * y);
            }
        }
        
        public double CylindricalRadialVelocity
        {
            get
            {
                double R = CylindricalRadius;
                return 1.0 / R * Math.Sqrt(vx*vx * x*x + vy*vy * y*y);
            }
        }
        
        public double CylindricalTangentialVelocity // ?
        {
            get
            {
                double p = Phi;
                return -1.0 * vx * Math.Sin(p) + vy * Math.Cos(p);
            }
        }
        
        public double CylindricalZVelocity
        {
            get
            {
                return vz;
            }
        }
        
        public double EllipsoidalRadius(double q, double s)
        {
            return Math.Sqrt(x * x + y * y / q / q + z * z / s / s);
        }
                
        public override string ToString()
        {
            return String.Format("{0:f6} {1:f6} {2:f6} {3:f6} {4:f6} {5:f6} {6:f6}",
                mass, X, Y, Z, Vx, Vy, Vz);
        }
        
        public int CompareTo(object o)
        {
            double r_this = RadiusSquared;
            double r_o = ((Particle)o).RadiusSquared;
            
            if (r_this > r_o)
                return 1;
            else if (r_this < r_o)
                return -1;
            else
                return 0;
        }
            
    }
}
